/* 
 * File:   bandpassfilter.hpp
 * Author: gian
 *
 * Created on May 7, 2012, 8:11 AM
 */

#ifndef BANDPASSFILTER_HPP
#define	BANDPASSFILTER_HPP

#include <gsl/gsl_complex_math.h>
#include "vector.hpp"

typedef gsl_complex cdouble;
typedef gsl_complex_float cfloat;
typedef double * doubleP;
typedef cdouble * cdoubleP;

class bandpassfilter {
public:    
  bandpassfilter(int N, double A, double fs, double fstop, double f0);
  virtual ~bandpassfilter();
  gsl_complex_float operator()(float x);
  void apply(cvectorf *out, vectorf *in);    
    
private:

  void Initialize_filter(int N);
  void Calculate_Chebichev(double A, double fs, double fstop, double f0);
  void Calculate_Lut(double f0);
  int Cheby2(int N, double A, double os, cdouble *zeros, cdouble *poles);
  void zero_pole_to_coefficients(cdouble zero, cdouble pole, double *a, double *b);
  void zero_pole_to_coefficients1(cdouble zero, cdouble pole, double *a, double  *b);
  cdouble double_order_filt(cdouble x, double *a, double *b, cdouble *w);
  cdouble first_order_filt(cdouble x, double *a, double *b, cdouble *w);
    
  int MN; /*Number of zeros/poles */
  int Mnb2; /* Number of second order blocks*/
  int Mnb1; /* Number of first order blocks*/
  double MA; /*Stopband attenuation */
  double Mfs; /*Sampling frequency */
  double Mfstop; /*Stopband edge */
  double Mos; /*Stopband edge in normalized units */
  double Mf0; /*Centra frequency of filter (rounded)*/
  gsl_complex *Mzeros; /*Zeros */
  gsl_complex *Mpoles; /*Poles */
  double *Ma1; /*First order filter coefficients */
  double *Mb1; /*First order filter coefficients*/
  double **Ma2; /*Second order filter coefficients */
  double **Mb2; /*First order filter coefficients */
  gsl_complex *Mw1; /*Internal register of first order block */
  gsl_complex **Mw2; /* Internal registers of second order blocks*/
  int MLutLen; /*lookup table length*/
  gsl_complex *MLut;
  long int Mindex;    
    
    
};

#endif	/* BANDPASSFILTER_HPP */

